Resistivity scaling and critical dynamics of fully frustrated Josephson-junction arrays 

with on-site dissipation 
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We study the scaling behavior and critical dynamics of the resistive transition in Josephson- 
junction arrays, at / = 1/2 flux quantum per plaquette, by numerical simulation of an on-site 
dissipation model for the dynamics. The results are compared with recent simulations using the 
resistively-shunted-junction model. For both models, we find that the resistivity scaling and critical 
dynamics of the phases are well described by the same critical temperature as for the chiral (vortex- 
lattice) transition, with a power-law divergent correlation length. The behavior is consistent with 
the single transition scenario, where phase and chiral variables order at the same temperature, but 
with different dynamic exponents z for phase coherence and chiral order. 
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I. INTRODUCTION 



There has been considerable interest, both ex- 
perimentally and theoretically, in phase transi- 
tions of two-dimensional Josephson junction arrays 

r j J ^) 1.2.3.4.5.6.7.8.9.10.11.12.13.14.15.16.17.18.19.20.21.22.23.24.25.26.1 ; 

Such arrays can be artificially fabricated as a lattice 
of superconducting grains connected by an insulator or 
a normal met&)l&2^£& and are also closely related to 
superconducting wire networksi^ Experimentally, the 
resistive transition has been the one most extensively 
studied(2i2iiiSiL& while theoretically several studies of XY 
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which describe the ideal JJA, have been done. A signif- 
icant understanding of these systems has already been 
achieved by comparing the results of experiments with 
the theoretically predicted equilibrium critical behavior, 
with and without an applied magnetic field. However, 
to a large extent, dynamical critical behavior remains 
much less understood, particularly in the presence of a 
magnetic field, where frustration effects may introduce 
additional elementary excitations relevant for the static 
and dynamic critical behavior. It is well known that 
while static critical phenomena depend on the spatial 
dimensionality as well as on the symmetry of the order 
parameter, the dynamic universality class of the phase 
transition will depend upon additional properties which 
do not affect the statics as, for example, conservation 
laws for the order parameter^ Thus, testing the univer- 
sality hypothesis of dynamical critical behavior requires 
the study of specific dynamical models. For JJA, the 
physically relevant dynamical model for the phase 
dynamics can not be unambiguously identified^iiSiiiiiS 
and should depend on the particular coupling mechanism 
between the superconducting elements of the array. It 
is expected that, at least for an array of ideal tunnel 



junctions, the Resistively-Shunted- Junction (RSJ) model 
of current flow between superconducting grains would 
be a more physical representation of the system^ 
This model assumes that energy dissipation occurs only 
through the junctions and imposes current conservation 
at each site of the array. On the other hand, for wire 
z**jia^i&.aa*^3fuxywMuu. effect junctions, local 

dissipation at the sites of the array should also be 
allowed leading to a model with on-site dissipation for 
the dynamics. 

In experimental investigations of JJA,, the resistive 
transition is usually identified from the behavior of the 
current- voltage (IV) characteristics near the critical tem- 
'^r ! ii' 3 i: 35 .' 36r ??j' 38 '??" 4( t*-nt correlation length determines 
both the linear and nonlinear resistivity sufficiently close 
to the transition. To interpret the data and determine 
the underlying equilibrium transition, the scaling the- 
ory of Fisher et idM- has been widely used. For JJA at 
zero magnetic field, which is isomorphic to the standard 
XY model, the resistive transition is in the Kosterlitz- 
Thoulcss (KT) universality class^i^i^ where the corre- 
lation length diverges exponentially near the critical tem- 
perature. Studies of the critical dynamical behavior using 
Monte Carlo (MC) dynamics^ and RSJ or on-site dissi- 
pation dynamics jiLiS find a behavior consistent with the 
dynamical theory of the KT transition. The exponent of 
the current- voltage relation, V ~ I a , at the transition, 
assuming the universal value a — z + 1 = 3, corresponds 
to a dynamic exponent z = 2 in the resistivity scaling 
theory^ 

However, in frustrated Josephson-junction arrays 
(FJJA), corresponding to / = 1/2 flux quantum per pla- 
quette, besides the phase variables, the vortex-lattice in- 
duced by the external field introduces an additional dis- 
crete (Ising-like) order parameter, the chiralityji^ which 
measures the direction of local current circulation in the 
array. The ground state consists of a pinned commensu- 
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rate vortex-lattice corresponding to an antiferromagnetic 
arrangement of chiralities and vortex-lattice melting cor- 
responds to the chiral order-disorder transition. As a 
consequence, two distinct scenarios for the occurrence of 
phase transitions as a function of temperature have been 
proposed by Teitel and Jayaprakash (Ref. flfilfb'l'l: sepa- 
rated chiral and phase-coherence transitions or a single 
transition where phase and chirality order at the same 
temperature. In the former scenario, the phase transi- 
tions should be in the KT and Ising universality classes, 
respectively, while in the later scenario the critical be- 
havior should be a superposition of KT and Ising critical 
behavior at the same critical temperature, if the cou- 
pling between phase and chiral variables are irrelevant at 
criticality (decoupled single transition), otherwise new 
critical behavior (coupled single transition) should oc- 
cur with phase coherence and chiral order showing crit- 
ical behavior different from the KT and Ising univer- 
sality classes. These possible scenarios are supported, 
for example, by renormalization-group studies based on 
the Ginzburg-Landau expansion of the frustrated XY 
model (FXY)ii which also shows that the universality 
class of these transitions can be described by the XY- 
Ising models. It appears that the current predominant 
point of view is that the separated transition scenario is 
realized with a KT transition occurring below the chi- 
ral transition. Recently, this scenario has received fur- 
ther support from appealing arguments by Korshunov 28 , 
based on chiral domain wall fluctuations and vortex un- 
binding, which provides a mechanism for the separation 
of the two transitions in this order. Also, there are signif- 
icant numerical evidences from equilibrium calculations 
which favor this scenario. However, the coupled single 
transition scenario has also received some support from 
different calculations of the chiral critical exponents and 
central charge from finite-size scaling which show results 
different from the pure Ising values, but several of these 
studies do not verify if the transition temperature for 
phase-coherence coincides with the chiral transition tem- 
perature. On the other hand, the numerical evidence for 
the separated transitions scenario finds that the phase- 
coherence transition 2 ^ is consistent with KT behavior but 
the critical exponents found for the chiral transition by 
finite-size scaling do not show the expected pure Ising be- 
havior according to Ref. However, it has been found 
by Olsson^i that the value of the thermal critical expo- 
nent is consistent with pure Ising value depending on the 
temperature region in which a fit is made. Therefore, 
the deviations of the exponents from pure Ising values 
can not be regarded as an unambiguous evidence for non- 
Ising critical behavior. The separated transition scenario 
also relies on the assumption that the phase-coherence 
transition is pure KT and therefore uses some of the 
predicted behavior from the KT theory, like the helic- 
ity modulus jump or exponentially divergent correlation 
length, to locate this critical temperature. If the helicity 
modulus jump is actually larger than the universal value 
then the procedure of locating the critical temperature 



from the jump 2 * can only overestimate the critical tem- 
perature. Although this assumption is consistent with 
a phase-coherence critical temperature below the chiral 
transition, such procedure could result in an underesti- 
mate of the phase-coherence critical temperature if the 
transitions coincide or the chiral transition occurs be- 
low the phase-coherence transition. In fact, it has been 
shown that if one enlarges the parameter space of the 
FXY models by considering a model where every other 
column in the square lattice has coupling constants which 
differ from the others by a constant ratio /?, the chiral 
transition occurs below the phase-coherence transition 30 
if p is sufficiently different from 1. It is then found that 
there is a singular contribution to the temperature depen- 
dence of the helicity modulus near the chiral transition 29 
determined by the chiral critical exponents. For the stan- 
dard FXY model, obtained when p — > 1, such singular 
contribution will remain if the transition is single and 
therefore it can affect the helicity modulus behavior near 
the transition. 

In any case, independently of the scenario interpreta- 
tion, several numerical calculations using quite different 
methods agreeSSiS&SliSii 3 ^ with the earlier estimate of the 
chiral transition temperature^ at T c h — 0.455 within 
a 0.8% errorbar. On the other hand, for the phase- 
coherence transition, it is clear that it would be more 
satisfactorily if it could be determined by methods which 
do no rely on assumptions of KT behavior. 

These different phase-transition scenarios have impor- 
tant consequences for the resistive behavior of the FJJA. 
Since the resistive transition corresponds to the onset of 
phase coherence, they imply quite distinct behavior. In 
the separated transition scenario or single but decoupled 
scenario, the resistive behavior should be described by 
the KT universality class. On the other hand, in the sin- 
gle coupled scenario, where the critical dynamics involve 
strongly coupled phase and chiral variables, the resistive 
behavior should be significantly different. In principle, 
such behavior can be detected experimentally. 

Measurements of current-voltage curves in FJJA were 
fitted assuming pure KT behavior but either an un- 
expectedly low value of the transition temperature (com- 
pared with theoretical expectations) was obtained in one 
case^ or the IV exponent at the transition was a < 3 in 
the other case 3 -. More recently, the current-voltage curves 
in JJA^ and in superconducting networks^ were found 
to be better described by a power-law correlation length. 
However, very different values of the critical exponents 
z, v were obtained in each case. 

Earlier numerical studies of the IV characteristics 
for FJJA, obtained with RSJ dynamics 3 ^ or MC 
dynamics^ were performed for small system sizes (L < 
16). In particular, the studies with RSJ dynamics 
used free boundary conditions to impose a driving cur- 
rent. This leads to significant additional dissipation due 
to boundary effects^ specially in small system sizes. 
Other works have studied the short-time dynamics of 
chiralityj 3 ^ and the non-equilibrium transitions at large 
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currents mL 

Recently^ we have studied the critical dynamics and 
resistivity scaling in FJJA by numerical simulation of the 
RSJ dynamics with periodic (fluctuating twist) boundary 
conditions including much large systems sizes. It was 
found that the current-voltage scaling is consistent with 
the single-transition scenario. The scaling behavior is 
well described by a resistive transition occurring at a crit- 
ical temperature corresponding to the chiral transition, 
with a power-law divergent correlation length, but with 
two different dynamic exponents, z p h ~ 1 and z c h ~ 2, 
for phase and chiral variables, respectively. This result 
implies that, at the transition, the exponent of the IV 
power-law, V ~ I a , is a = z P h + 1 ~ 2 rather than a = 3 
as for the unfrustrated case. In view of the possible de- 
pendence of the dynamic behavior on the particular RSJ 
dynamics used in these simulations, it should also be of 
interest to study the resistive behavior with an on-site 
dissipation model for the dynamics. Results for this dy- 
namical model should be particularly relevant for frus- 
trated wire networks^ or proximity-effect junctions. 

In this work we study the resistivity scaling and crit- 
ical dynamics of a frustrated Josephson-junction array, 
defined on square lattice, at / = 1/2 flux quantum per 
plaquette, by numerical simulations of an on-site dis- 
sipation model for the array dynamics. Using a dy- 
namic scaling analysis, we find that the resistivity be- 
havior and critical dynamics are well described by the 
critical temperature corresponding to the chiral (vortex- 
lattice) transition with a correlation length that diverges 
as a power law. Two dynamic exponents, z p h ~ 1.5 and 
z c h ~ 2.5, are found for phase-coherence and chiral or- 
der, respectively. Consequently, at the transition, the 
exponent of the current- voltage power-law, V ~ I a , is 
& = z.ph + 1 ks 2.5 rather than a = 3 as for the un- 
frustrated case. This is the same behavior we have found 
recently for the RSJ models but with different values for 
the dynamic exponents (z p h ~ 0.9(1) and z c h ~ 2.1). In- 
cluding on-site dissipation in the dynamical model could 
be a more realistic description of wire networks than the 
RSJ model. Indeed, resistivity scaling of experimental 
data on wire networks 8 find z ~ 2, which is consistent 
with our estimate z p h within the experimental errors, and 
also shows that the resistivity scaling is well described 
by a power-law correlation length as found in our simu- 
lations. 



II. MODEL AND SIMULATION 

The hamiltonian of a square two-dimensional array un- 
der a magnetic field is given by 

H = -Ej cos(0 r+M -0 r - A r ,„) (1) 

where 8 r is the superconducting phase of the grain at 
site r = {n x a,n y a) with n x ,n y integers, and a the lat- 
tice constant, and // = x, y with x = (a, 0), y = (0,a), 



and Ej = Ioh/2e the Josephson energy. The magnetic 
field introduces frustration through the vector potential 
integral 

2?r r T+f * 

which satisfies 

= 2tt/, (3) 

with / = Ha 2 /$Q, where H is the applied magnetic field 
and <&o = h/2e the quantum of flux. The fully frustrated 
case corresponds to half quantum of flux per plaquette, 
/ = l/2. 

The simulations are performed with the same "fluctu- 
ating tw ist" boun dary conditions as used, for example, 
m Refs. Illl39l4fl This consists on considering periodic 
boundary conditions for the supercurrents in the \x di- 
rection while adding a fluctuating twist a M to the gauge 
invariant phase in the fi direction. In this case the gauge 
invariant phase difference is modified to 

0r,/i = 0r +t i - Or- A r ,„ + d p (4) 

For the vector potential we choose the Landau gauge 
A r x = —2irfriy 

A I<y = (5) 

In this gauge, the boundary condition for the phases in 
a system of size L x L is given by 

6{n x + L, n y ) = 9(n Xl n y ) 

9(n x ,n y + L) = 9(n x ,n y ) - 2nfLn x . (6) 

For / = 1/2 and L even, the second condition is irrele- 
vant, but not for general frustration /. In the presence 
of an external current J^ xt in the /x direction, one has 
to add the term —^L 2 I^ t a^ in the hamiltonian of Eq. 
JTJ, which couples the current with the global phase dif- 
ference per row, La^, introduced by the fluctuating twist. 
Therefore, the hamiltonian of a frustrated square array 
with fluctuating twist boundary conditions and an exter- 
nal current is 

H = -Ej cos((9 r+M -9 t - A T ^ + a^) 

-^'E^ (7) 

It 

We define the on-site dissipation dynamics by consider- 
ing the local Langevin equations for the fluctuating vari- 
ables 9 r and 

w = < 8 > 
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where To, T a are dissipation parameters, and the noise 
terms have zero average and correlations 

(Vr(t)T}r'(t')) = 2k B TTg6r,r'S(t-t r ) (10) 

(%(«)v(0) - 2k B TT a S^S(t-t') (11) 

The dissipation constant r Q should be proportional to 
L~ 2 in order to be an intensive quantity. A convenient 
choice is 



(in general it can be r Q = 3Tg/L 2 , here we choose (3—1 
to be consistent with Ref . Il2|) . 

Dimensionless quantities are used with time in units of 
t = 2e/hTgIo, currents in units of Iq, voltages in units 
of (h/2e) 2 TgI and temperature in units of TiI /2ekB- A 
total current / is imposed uniformly in the array in the 
y-direction with current density J = I/L, where L is the 
system size and the average electric field E is obtained 
from the voltage V across the system as E — V/L = 
(h/2e)(da y /dt), where a y L is the global phase difference 
or twist in the y-direction. With all this considerations, 
the dimensionless equations of motion are then 

^ = -A,-S r ,,+ Vr (t) (12) 

r 

where the supercurrent is defined by 

S rili = sm(0 r+fl -6 T - A Ttll + otft) , 
and the discrete divergence operator is defined as 

fj.=x,y 

Finally, the now dimensionless noise variables rj r {t) have 
correlations 

(Vr(t)Vr'(t')) = 2T8r,r'5(t-t) (14) 

2T 

(1,WV(0) = jj6^6(t-t). (15) 

The set of Eqs ©-© describe the dynamics of JJA 
with "on-site dissipation" in contrast to the RSJ dy- 
namics which only considers dissipation through the 
junctions 13 . The on-site dissipation dynamical model 
has been studied previouslySiii* 1 ^ . for the unfrustrated 
(/ = 0) case and compared with the RSJ dynamics. 
Their main difference is that while the on-site dynamics 
corresponds to a local damping the RSJ dynamics corre- 
sponds to a nonlocal damping 9 *^. A physical interpreta- 
tion of the on-site dynamics for JJA in terms of currents 
and voltages has also been discussed previousl y 9 ' 11 ' 12 ' 40 . 
Its main features are summarized in the following, (i) It 
takes into account normal current flow between each su- 
perconducting node and the substrate, which leads to 



a current leakage through a resistance to the ground 
Rq. (ii) It neglects the quasiparticle normal current of 
each junction, which is associated with a shunt resistante 
R s . This means taking R s — > oo, or actually assuming 
R s 3> Rq for the array. The assumptions (i) and (ii) lead 
to Eq. © for 6 T , which corresponds to the conservation 
of supercurrents at each node plus a leakage of normal 
current to the substrate. In this case we get 

T e = (2e/h) 2 R . 

However, if one considers Eq. JHJ alone for the calculation 
of current- voltage curves with open boundary conditions 
it is found that an applied external current leads to dis- 
sipation only at the boundaries were current is injected 
(extracted), since normal current will flow directly from 
the first (last) row of junctions to the substrate through 
Strictly periodic boundary conditions are not pos- 
sible to be implemented in a consistent way. (iii) In order 
to correctly model current-voltage curves and to be able 
to implement fluctuating twist periodic boundary condi- 
tions, one has to addii* 4 ^ a global normal current channel 
in parallel to the whole array, with a "global resistance" 
i? g iobai, such that in the normal state the total resistance 
of the array will be given by -ffgiobai- Then total con- 
servation of current leads to Eq. @ which represents a 
parallel circuit of the average supercurrent in the array 
and the global normal current. In the approach of Refs. 
Illlll2ll40l R„]„h*] = Rq is assumed, and therefore this leads 
to the choice T a = (2e/h) 2 R /L 2 . 

The above assumptions (i)-(iii) for the on-site dynam- 
ics give a consistent interpretation of calculations of the 
current-voltage response and phase dynamics, but corre- 
spond to a model system rather than to a particular JJA 
available experimentally. In realistic JJA the normal cur- 
rents in the junctions can not be neglected since usually 
R s <C i?o, and therefore the RSJ model can be a good 
representation of the JJA. A possible realization of the 
dynamics of Eqs. ©-© could be achieved experimen- 
tally if one adds on purpose a resistor in parallel to the 
whole array such that it has a resistance i? g i bai •C R s , 
in which case normal currents will mainly go through 
^global and reduce the weight of the normal quasiparticle 
currents of the junctions 4 ^. 

A good candidate for the on-site dynamics is a super- 
conducting wire network. In this case one has to take into 
account the dynamics of the complex order parameter 
which is given by the time dependent Ginzburg-Landau 
equation (TDGL) coupled to the electromagnetic field 
equations 4 ^ 4 ^ 4 ^. There are two dissipative mechanisms 
in this case: (i) via the normal state resistivity, since the 
total current is the sum of the supercurrent and the nor- 
mal current in each wire of the network (this is the equiv- 
alent of the shunt resistance of the RSJ model in a JJA), 
and (ii) via the relaxation of the complex order param- 
eter in the TDGL equations, which is local in nature 4 ^ 
and where its dominant contribution is determined by 
D, the normal state diffusion constant. After writing 
the TDGL equations in a discrete lattice, and neglecting 
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the fluctuations of the amplitude of the order parameter 
(London limit), one obtains^ that the on-site part in the 
dynamics of the phase is provided by D, and would cor- 
respond to Eq. © with Tg = 16tt 3 DA 2 /($%aS), where 
a is the network lattice constant and S the section of 
the wires. Therefore, there is no need to invoke a "leak- 
age of normal current to the ground" in this case. The 
full dynamics of the superconducting wire network is a 
mixture of both the "on-site" dynamics and the "RSJ" 
dynamics. However, in the presence of an on-site con- 
tribution, the resulting rate of change of the phases at 
different sites, like eq. (8), does not have a logarithmic 
nonlocal dependence at large separations as in the pure 
RSJ models 

In any case, in the present work, we will take the purely 
on-site dynamical equations as a model dynamics that 
corresponds to a limit of the general dynamics of a JJA 
or a superconducting wire network where only local dis- 
sipation is taken into account. The opposite limit for the 
dynamics is the pure RSJ model that we have analyzed 
in Ref. S 

We integrate the dynamical equations with a second 
order Runge-Kutta-Hclfand-Grccnside method with time 
step At — 0.01 — 0.07t, averaging over, typically, 10 6 time 
steps after using 5 x 10 5 time steps for equilibration. The 
results were averaged over 5 — 10 different initial configu- 
rations of the phases and system sizes ranging from L = 8 
to L = 180 were considered. 



The scaling form of Eq. ^]does not take into account 
finite-size effects and so it is valid only in a range of 
temperature T and current densities J where such effects 
are not dominant. Finite-size effects are very important 
sufficiently close to the transition when the correlation 
length £ reaches the system size L. In particular, at T c , 
the correlation length £ will be cut off by the system size 
L in any finite system. From Eq. (JTSJ, the nonlinear 
resistivity at T c should then satisfy the scaling form 



(17) 



It follows from Eq. (|17|) that the linear resistance Rl = 
limj_>.o E/J should decrease as a power- law of the system 
size, 



Rl cx L 



(18) 



right at T c . This behavior is independent of the form of 
the correlation length divergence. The linear resistance 
can be obtained from the Kubo formula of equilibrium 
voltage fluctuations as 



Rl = 



1 

2T 



dt(V(t)V(0)). 



(19) 



without an imposing driving current. Rl can also be 
determined more accurately from the long-time fluctu- 
ations of the total phase difference across the system 
A6{t) = La„ asiS-ii 



III. DYNAMIC SCALING THEORY 

Near a second-order phase transition, the diverging 
correlation length £ leads to critical slowing down charac- 
terized by relaxation times r that also diverge approach- 
ing the transition temperature. The dynamic scaling 
hypothesis^ asserts that measurable quantities should 
scale with the diverging correlation length £ and the re- 
laxation time t cx £ z , near the transition temperature, 
where z is the dynamical critical exponent. A general 
dynamic scaling theory for the resistivity behavior near a 
superconducting transition has been provided by Fisher, 
Fisher, and Huse4£ According to this scaling theory, the 
nonlinear resistivity E/J should satisfy the scaling form 



(16) 



in two dimensions, where the + and — correspond to 
the behavior above and below the transition, respec- 
tively. For a transition in the KT universality class, 
the correlation length should diverge exponentially as 
£ cx exp(6/|T/T c - l\ 1/2 ), above T c . Otherwise, for a 
usual continuous transition, a power-law behavior is ex- 
pected, £ cx \T/T C — 1| _1/ , with an exponent v to be 
determined. Thus, a scaling plot according to Eq. i|16|) 
can be used to verify the dynamic scaling hypothesis and 
the assumption of an underlying equilibrium transition. 



RL = —({A6(t)-A6{0)Y) 



(20) 



valid for sufficiently long times t. 

The critical dynamics leading to the resistivity scaling 
described above can also be studied by the behavior of 
time correlation functions. For the frustrated JJA, there 
are two different types of time correlations of particular 
interest, the time correlation for chiralities C c h(t) and 
phase variables C p h(t). We shall use normalized time 
correlation functions defined as 



C{t) 



(A(t)A(0)) {Af 
(A 2 ) (A) 2 



(21) 



For the phase variables, A = S = ^2 { Si, where s = 
(cos(0), sin(#))) and for the chiral variables A = \ — 
J2<ij>(@i ~ @j ~ Aij)/2w, where the summation is taken 
over the elementary plaquette of the lattice and the 
gauge-invariant phase difference is restricted to the in- 
terval [— 7r,7r]. The relaxation time r can be obtained 
from the exponential decay C(t) cx exp(— t/r) at suffi- 
ciently long times. In general, the time dependence of 
C(t) can be expressed as a series of exponential terms 
with the largest decay time corresponding to the critical 
relaxation time of the long time dynamics^. From dy- 
namic finite-size scaling, the relaxation time should scale 
at T c as r cx L z , from which the z can be estimated from 
the slope in a log log plot. An alternative procedure to 
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estimate z from equilibrium dynamics is to explore the 
expected finite-size behavior of the time correlation func- 
tions at long times. Since at T c the relaxation time scales 
as t (x L z , the time correlation function for different sys- 
tem sizes can be cast into a scaling form in terms of the 
dimensionless ratio t/L z as 

C(L,t) = C(t/L z ) (22) 

where C(x) is a scaling function. However, this assumes 
a simple scaling form for the time correlation functions 
and is only valid for sufficiently long times when a single 
exponential term describes the relaxation behavior. 



IV. RESULTS AND SCALING ANALYSIS 

Fig. n shows the temperature dependence of the non- 
linear resistivity E/J for the largest systems size L = 180 
near the chiral transition temperature T c h, estimated 
previously from equilibrium Monte Carlo simulation, 18 
T c h = 0.455. Qualitatively, the linear resistance Rl = 
limj^o E/J, tends to a finite value at high temperatures 
but extrapolates to very low values at lower tempera- 
tures, consistent with the existence of a resistive transi- 
tion in the range 0.45 < T c < 0.46. In the double tran- 
sition scenario, where the phase-coherence transition is 
expected to be in the KT universality class, the estimate 
of the proposed KT critical temperature is Tkt = 0.446, 
from Monte Carlo simulations^ which is close to T c h- 
However, as it is clear from the behavior at the lowest cur- 
rents in Fig. this estimate is below the resistive transi- 
tion since the resistivity curves for T = Tkt = 0.446 and 
T = 0.45 > Tkt tends to zero for J — > 0, indicating that 
the system is still in the superconducting phase for these 
temperatures. On the other hand, the resistivity curve 
for T — 0.46 > T c h clearly tends to finite resistivity for 
J — » 0. This shows that the resistivity transition occurs 
at T c h or at a temperature very close to T c h rather than 
at the proposed estimate of Tkt- 

Additional support for a resistivity transition at T c h = 
0.455 is provided by the behavior of the linear resistiv- 
ity Rl as a function of system size, shown in Fig. [21 
For T > 0.455, Rl extrapolates to a finite value con- 
sistent with the behavior of the nonlinear resistivity for 
J -> in Fig. [U On the other hand, for T < 0.455 it 
extrapolates to zero, indicating that the resistive tran- 
sition temperature is compatible with the estimate of 
T c h = 0.455. Since in this calculations Rl is obtained 
without any current bias, from the equilibrium dynam- 
ical fluctuations, according to Eq. ^5] this result also 
verify that the T c inferred from the behavior of the non- 
linear resistivity for the largest system size in Fig. ^ is 
not an artifact of finite current bias and in fact reflects 
the underlying equilibrium transition behavior. 

Although the resistivity behavior of Figs. Hand HI al- 
ready suggest that the resistive transition temperature 
coincides with T c h or it is much closer to this value than 



previous estimates, we now proceed, as in any study of 
critical phenomena, to obtain the asymptotic equilibrium 
critical behavior in the thermodynamic limit, L — » oo 
and J — > 0, from a scaling theory. A scaling plot accord- 
ing to Eq. ^]is shown in Fig. [3] for the largest system 
sizes, in the temperature range closest to T c h and small- 
est current densities, assuming the correlation length £ 
has a power-law divergence with T c — T c h and using v 
and z as adjustable parameters so that the best data col- 
lapse is obtained. This scaling plot shows that the two 
largest system sizes L = 128 and L = 180 give the same 
data collapse and so finite size effects neglected in the 
scaling form of Eq. 1161 are not dominant for the range 
of temperatures and current densities shown in the plot. 
Similar scaling analysis assuming a KT correlation length 
and fixing T c at the estimate of Tkt does not result in 
a good data collapse. The same behavior was found us- 
ing the RSJ dynamics^ 8 - From this scaling analysis, we 
estimate v = 0.9(1) and the dynamical critical exponent 
z = 1.3(3). The static exponent v is consistent with es- 
timates of the chiral transition from equilibrium Monte 
Carlo simulations^ but the accuracy is not sufficient to 
rule out the value v = 1 expected for the standard Ising 
transition. Our estimate of z is smaller than the one ob- 
tained previously for the frustrated XY model with MC 
vortex dynamics^ where z ~ 2 was found. However, 
such MC simulation corresponds to a different dynamics 
and also only very small system sizes (with L = 8 — 14) 
were analyzed. We now take into account finite-size ef- 
fects explicitly by studying the scaling behavior of the 
linear resistivity Rl near T c in Fig. [21 At T c , the lin- 
ear resistivity should scale with system size according to 
Eq. ^| Near T c , it should also depend on temperature 
through the dimensionless variable £/£. If the correla- 
tion length diverges as a power law then it should satisfy 
the finite-size scaling form 

R L L' = }{(T/T C - 1)L^) (23) 

In fact, as shown in Fig. 0] the linear resistivity data 
satisfy the scaling form with T c = T c h and a value z = 
1.5(2) consistent with the estimate from the nonlinear 
resistivity scaling. 

The above scaling analysis for the nonlinear resistiv- 
ity at large system sizes and linear resistivity at smaller 
system sizes already confirm that the resistive transi- 
tion temperature T c is very close to T c h, with a dy- 
namic exponent z < 2. However, in the absence of a 
completely satisfactorily determination of T c from static 
critical behavior^2*2iiS4 from now on, we will assume 
T c = T c h and explore to which extent this give us consis- 
tent results for the dynamical critical behavior, including 
finite-size effects. Another reason to assume the value of 
T c obtained from equilibrium simulations rather than es- 
timating from the dynamic scaling itself is that, in gen- 
eral, the most reliable way of studying critical dynam- 
ics and determine the dynamic exponent z is to use the 
known value of T c . This is true not only for models where 
T c is known exactly as for the two-dimensional Ising 
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model^i but also for models where T c is only known by 
numerical simulations as for the three-dimensional Ising 
models 

An alternative estimate of z can be obtained from the 
nonlinear resistivity by studying the expected size de- 
pendence at T c . As shown in Fig. [SJ a finite size scaling 
according to Eq. (|17|l gives the same dynamic exponent 
z = 1.4(3), within the estimated error bar. The same be- 
havior was also observed using the RSJ dynamics^ but 
with a smaller value of z. Equilibrium calculations of the 
linear resistance Rl at T c h also give a consistent estimate. 
Fig. [SJshows the finite size behavior of Rl obtained from 
Eq. (|18|l . A power-law fit gives z = 1.41(5) which agrees 
with the other estimates and suggests therefore that the 
value of z corresponds to the underlying equilibrium dy- 
namical behavior. To show the reliability of this method, 
it is also included in Fig. [B]the behavior for the unfrus- 
trated case, / = 0. In this case the resistive transition 
is in the KT universality class and a dynamical expo- 
nent z = 2 is expected, independent of the dynamics. 
Indeed, for / = 0, the same power-law fit at the criti- 
cal temperature T c = 0.887 estimated from Monte Carlo 
simulations^ gives z — 1.96(5), in good agreement with 
previous resistivity calculations^ for / = using smaller 
system sizes up to L = 16. 

It should be noted that our above estimate of the dy- 
namic exponent z is obtained by requiring that T Cl z and 
v satisfy at the same time the finite-size scaling forms of 
Eqs. i|17[l . 1|18|) and i|23|) . including small system sizes, 
as well as the scaling form of Eq. I|16fl f° r the largest 
system sizes. Using only Eq. (|16f) can lead to inaccu- 
rate estimates of z as shown recently in Ref. |53j for the 
unfrustrated case. 

To further verify that the estimate of z obtained from 
the resistivity scaling does in fact reflect critical phase 
fluctuations near the transition rather than just criti- 
cal fluctuations for the chiral order parameter, we have 
also performed equilibrium calculations of the phase au- 
tocorrelation functions C p h (t) for the phase variables and 
C c h(t) for the chirality variables. Figs. [7] and [S] show 
the finite-size behavior of the time correlations functions 
evaluated at the critical temperature T c h- If this temper- 
ature corresponds to the critical point for phase coher- 
ence and vortex-lattice disorder then the relaxation times 
for both phase and chirality variables should diverge with 
the system size as t a i 2 . The relaxation time T p h and 
T c h can be obtained from the exponential decay of C(t) 
at sufficiently long times. We take into account possi- 
ble contributions from short time behavior by fitting the 
time dependence of C p h(t) and C c h(t) to a sum of two 
exponentials and extract r from the largest decay time. 
Fig. [5] shows the finite-size behavior of the relaxation 
time at T c h for the phases and chiralities. From a power- 
law fit we obtain z p h — 1-8(1) from the phase relaxation 
time T p h which is indeed consistent, within the estimated 
error bar, with the value of z obtained from the resis- 
tivity scaling discussed above. The estimate from the 
chiral relaxation time in Fig. [5] is significantly different, 



Zch = 2.5(2). For an alternative estimate of z we have 
also used the scaling of the correlation function itself. 
The time correlation functions should satisfy the scaling 
behavior of Eq. [221 As shown in Figs. HUlandlTTl C p h(t) 
and C c h {t) indeed satisfy the expected finite size behavior 
at the critical temperature providing additional estimates 
of the dynamic exponents z p h = 1-9(2) and z c h = 2.6(2) 
which are consistent, within the estimated error bar, with 
the values obtained from the relaxation time scaling. Fi- 
nally, above T c , the relaxation time should depend both 
on system size and temperature. If the correlation length 
diverges as a power law then r p h and r c h should satisfy 
the finite-size scaling form 

tL- z = f((T/T c - 1)LV") (24) 

In fact, the data collapse in Figs. 1121 and IT31 show that 
this scaling form is satisfied with T c = T c h and the val- 
ues of Tph and r c h which are consistent with the above 
estimates. 



V. DISCUSSION 

Recently, Holzer et alm& showed that for the unfrus- 
trated case, / = 0, the scaling behavior in Eq. i(Po|) 
considered alone, i.e., without taking into account finite- 
size effects, yields incorrect values for the dynamic ex- 
ponent z, using approximate analytical expressions for 
the IV characteristics available in the literature 1 . We 
should emphasize that our approach for the resistivity 
scaling analysis described in the previous Section is quite 
different. Our estimate of the dynamic exponent z is ob- 
tained by requiring that T c , z and v satisfy at the same 
time the finite-size scaling forms of Eqs. (|17fl . I|18|) and 
(|23|l . including small system sizes, as well as the scaling 
form of Eq. (|16|) for the largest system sizes. It should 
also be considered that the possibility of an equilibrium 
KT transition for / = 1/2 within the separated transi- 
tions scenario does not imply that the dynamics would 
be the same as the KT dynamics and therefore for the 
frustrated case considered here there is no reliable analyt- 
ical expressions available for the IV characteristics. The 
dynamics for f — 1/2 will be different because besides 
vortex excitations, chiral domain walls also contribute to 
the nonlinear resistivity as shown in Ref. |33[ Moreover, 
it has already been shown for the / — case that, when 
finite-size scaling is taking into account in the resistivity 
scaling theory of Fisher et alm&, as we also do in our ap- 
proach, the correct dynamic exponent z = 2 is obtained 
for the KT transition, as shown for example in Ref. H2T 
This is also verified in the scaling analysis of our data 
as shown in Fig. [S] where we find a dynamic exponent 
consistent with z = 2 for / = 0, as expected. 

The distinct values obtained for z p h and z c h with the 
on-site dissipation model deserve some considerations. 
Similar behavior was also found by us using the RSJ 
dynamics^ The final results for both models, obtained 
from the resistivity scaling and time-correlation function 
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scaling analysis, are summarized in Table 1. Although for 
the on-site model, the two methods give results for z p h 
which differ beyond the estimated errorbar, the values are 
significant below the value obtained for z c h- Naively, if 
the two transitions happen at the same temperature, one 
would expect that the same dynamic exponent should 
hold for the phase and chiral relaxation times. However, 
we should mention that different dynamic exponents for 
coupled order parameters have already been found previ- 
ously at multicritical points in magnetic systems^ 4 . This 
suggests that a possible explanation for two dynamic ex- 
ponents at the transition of the FJJA may rely on the ex- 
istence of a multicritical point in the phase diagram of the 
relevant effective Ginzburg-Landau free energy describ- 
ing the transition. A multicritical point is known to occur 
in the coupled XY-Ising modeli^ which should describe 
the static critical behavior of the FJJA and this could be 
a useful framework for investigations of the dynamical 
universality class of FJJA. In the context of supercon- 
ducting systems, different dynamic exponents for the re- 
sistivity and chirality have also recently been found in the 
resistive transition of disordered superconductors^ de- 
scribed by the three-dimensional XY spin glass models. 
Just as in the case of the frustrated JJA, the phase 
transition in the XY spin glass results from the com- 
petition of a chiral order parameter and phase variables. 
Although earlier work for this problem concluded for a 
spin-chirality decoupling picture^, more recent numeri- 
cal work have provided strong evidenced that there is a 
single transition at which both phase variables and chi- 
ralities order. 

Although the single transition scenario provides a con- 
sistent interpretation of our data, it is worth emphasiz- 
ing that the alternative separated transitions scenario^ 
can not be ruled out. We believe, there are two possi- 
ble explanations for some of our findings within the later 
scenario, as discussed below. 

It is possible that the KT transition is actually much 
closer to T c h than estimated previously and so the tran- 
sitions can not be resolved within the accuracy of our 
data. Our analysis of the resistivity behavior suggests 
that in this case it should occur above T c ~ 0.452. This 
value is already close or within the range of the error- 
bars reported for the chiral transition critical tempera- 
ture obtained, for example, by Monte Carlo simulations 
which gives T ch = 0.455(2) ( Ref. GJ) or 0.454(2) ( Ref. 
|23|) . It should be noted however that this only consid- 
ers the critical temperatures alone and not the critical 
behavior. In the alternative decoupled single transition 
scenario, the critical behavior should be described by a 
superposition of a pure KT and pure Ising transitions 
at the same critical temperature. However, this is also 
not consistent with our results. Nevertheless, even if the 
transitions are so close that their critical temperatures 
can not be resolved by any method, in principle, it could 
still be possible to distinguish these scenarios due to the 
mechanism discussed in Ref. [28] or due to the effects of 
different corrections to scaling. 



A second possibility is that the dynamic scaling the- 
ory of Fisher et alm-± in its original form in Eq. ^j] is 
not valid for the present case and should be enlarged to 
include the interplay of two divergent length scales at 
nearby temperatures 24 which can lead to crossover ef- 
fects at small length and time scales. In fact, the under- 
lying assumption in the resistivity scaling theory is that 
there is a single divergent length scale, corresponding to 
the leading divergent contribution to finite correlation 
lengths, when approaching the critical temperature of the 
resistive transition. This would certainly be valid within 
the coupled single transition scenario, which is consistent 
with our conclusions since in that case phase coherence 
and chiral order develop at the same critical temperature, 
with strongly coupled order parameters, and the equilib- 
rium critical behavior should be described by a single 
divergent length scale. Above the transition, in the dis- 
ordered phase, the chiral and phase correlation lengths 
diverge when approaching T c with a common leading di- 
vergent contribution. Below the transition, where there 
is chiral order and a Gaussian fixed line is expected for 
the phase variables, the chiral correlation length diverges 
when approaching T c with the same leading divergent 
contribution while the phase correlation length remains 
infinite since the Gaussian fixed line corresponds to the 
absence of a length scale. However, if phase coherence 
and chiral order develop at different temperatures then 
the resistivity scaling can only hold sufficiently close to 
the phase coherence transition otherwise the scaling form 
of Eq. 1 161 should be enlarged to include the divergent chi- 
ral correlation length in addition to the phase correlation 
length. This would lead to a scaling function g±(x, y) in 
Eq. depending on 2 scaling variables x — J£,kt /T 
and y = (, c h/£,KT, which makes the scaling analysis of 
the data very complicated specially when taking into ac- 
count finite-size effects. This could explain, for example, 
why a good scaling collapse like Fig. [SJis not obtained by 
assuming a resistive transition at T c = Tkt, estimated 
by previous works. However, it would remain unclear to 
us in this case why the linear and nonlinear resistivity 
scaling as well as the critical dynamics including differ- 
ent temperatures and system sizes are so well described 
by a resistive transition at T c = T c h- 



VI. CONCLUSIONS 

We have studied the resistivity scaling and critical 
dynamics of a frustrated Josephson-junction array, at 
/ = 1/2 flux quantum per plaquette, by numerical simu- 
lations of an on-site dissipation model for the array dy- 
namics. Using a dynamic scaling analysis, we find that 
the resistivity behavior and critical dynamics are well 
described by the critical temperature corresponding to 
the chiral (vortex-lattice) transition with a correlation 
length that diverges as a power law. Two dynamic ex- 
ponents, z p h ~ 1.5 and z c h ~ 2.5, are found for phase- 
coherence and chiral order, respectively. Consequently, 
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at the transition, the exponent of the current-voltage 
power-law, V ~ I a , is a — z p h + 1 « 2.5 rather than 
a = 3 as for the unfrustrated case. The same behav- 
ior has been found recently for the resistively-shuntcd- 
j unction model22. but with different values for the dy- 
namic exponents ( z p h ~ 0.9(1) and z c h ~ 2.1). One 
implication of these results for transport experiments 
is that the usual method of locating the critical tem- 
perature from the value corresponding to a nonlinear 
IV exponent a = 3, may lead to a significant under- 
estimate. This is more severe for tunnel-junction ar- 
rays which should be better described by the resistively- 
shunted-junction model^ where we expect a ~ 2 at the 
resistive transition^ For wire networks the on-site dis- 
sipation model should be more appropriate. Indeed, re- 
sistivity scaling of experimental data on wire networks 8 
find z ~ 2, which is consistent with our estimate of z p h 



within errors. It also shows that the resistivity scaling 
is well described by a power-law correlation length as 
found in our simulations. Further detailed IV measure- 
ments combined with magnetic properties, which could 
in principle probe the chiral transition, are needed to test 
our results. 
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TABLE I: Dynamic exponents of the resistive chiral transi- 
tion at T c h using the on-site dissipation model (TDGL) and 
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and C correspond to results obtained from the resistivity scal- 
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FIG. 1: Nonlinear resistivity E/J as a function of tempera- 
ture for system size L — 180. 
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FIG. 2: Linear resistance Rl, obtained without current bias, 
as a function of temperature and system size. Lines are just 
guide to the eyes 
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FIG. 3: Scaling plot of the nonlinear resistivity data for the 
smallest current densities near T c = T c h = 0.455 with £ oc 
\T/T C — Open symbols correspond to L = 128 and filled 

ones to L — 180. 
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FIG. 4: Finite-size scaling plot of the linear resistance data 
near T c = T ch = 0.455. 
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FIG. 5: Scaling plot of the nonlinear resistivity E/J at T c = 
T c h — 0.455 for different system sizes L. 




FIG. 6: Linear resistance as a function of system size at the 
critical temperatures T c = T c h f° r / = 1/2 and T c — 0.887 for 
f — 0. Power-law fits give estimates of the dynamic exponent 
z. 
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FIG. 7: Time correlation function C c h{t) for the chiral vari- 
ables at T c h, for different system sizes. 
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FIG. 8: Time correlation function C p h(t) for the phase vari- 
ables at T ch , for different system sizes. 
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FIG. 9: Finite-size behavior of the phase and chiral relaxation 
times, T p h and t c h respectively, at the critical temperature 
T c = T c h. Power-law fits give estimates of the dynamical 
exponents z p h and z c h- 




FIG. 10: Finite-size scaling plot of the time correlation func- 
tion C ch (t). 
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FIG. 11: Finite-size scaling plot of the time correlation func- 
tion Cph(t). 
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FIG. 12: Finite-size scaling plot for the relaxation time t c k 
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FIG. 13: Finite-size scaling plot for the relaxation time t p h 



